Identifying Vibrations that Control Non-adiabatic Relaxation of Polaritons in Strongly Coupled Molecule–Cavity Systems

The strong light–matter coupling regime, in which excitations of materials hybridize with excitations of confined light modes into polaritons, holds great promise in various areas of science and technology. A key aspect for all applications of polaritonic chemistry is the relaxation into the lower polaritonic states. Polariton relaxation is speculated to involve two separate processes: vibrationally assisted scattering (VAS) and radiative pumping (RP), but the driving forces underlying these two mechanisms are not fully understood. To provide mechanistic insights, we performed multiscale molecular dynamics simulations of tetracene molecules strongly coupled to the confined light modes of an optical cavity. The results suggest that both mechanisms are driven by the same molecular vibrations that induce relaxation through nonadiabatic coupling between dark states and polaritonic states. Identifying these vibrational modes provides a rationale for enhanced relaxation into the lower polariton when the cavity detuning is resonant with specific vibrational transitions.

1 Multi-scale molecular dynamics simulation model 1

.1 Extension of the Tavis-Cummings Hamiltonian to molecular systems
The multi-scale Molecular Dynamics (MD) method for simulating molecules under strong lightmatter coupling was presented previously. [1][2][3] For completeness, we provide a concise overview here. The method is based on the Tavis-Cummings Hamiltonian of quantum optics, and models the interaction between N molecules and n mode confined light modes in a one-dimensional Fabry-Pérot cavity ( Figure S1): [3][4][5][6] Here,σ + j (σ − j ) is the operator that excites (de-excites) molecule j from the electronic ground (excited) state |S j 0 (R j ) (|S j 1 (R j ) ) to the electronic excited (ground) state |S j 1 (R j ) (|S j 0 (R j ) ); R j is the vector of the Cartesian coordinates of all atoms in molecule j, centered at z j ;â kz (â † kz ) is the annihilation (creation) operator of an excitation of a cavity mode mode with wave-vector k z ; f z (z j ) = e ikzz j is the function describing the form of the quantized electromagnetic (EM) field modes, here taken to be that of plane waves with in-plane momentum k z ; hν j (R j ) is the excitation energy of molecule j, defined as: with V mol S 0 (R j ) and V mol S 1 (R j ) the adiabatic potential energy surfaces of molecule j in the electronic ground (S 0 ) and excited (S 1 ) state, respectively. The last term in Equation 1 is the total potential energy of the system in the absolute ground state (i.e. with no excitations in neither the molecules nor the cavity modes), defined as the sum of the ground-state potential energies of all molecules in the cavity.
The V mol S 0 (R j ) and V mol S 1 (R j ) adiabatic potential energy surfaces can be modeled with ab initio, density functional theory (DFT), or hybrid Quantum Mechanics / Molecular Mechanics (QM/MM) methods. While these methods provide access to both ground and excited state potential energy surfaces as well as transition dipole moments for complex systems with many molecular degrees of freedom, the accuracy depends on the level of theory in combination with the size of the atomic basis set. High accuracy results therefore require very large computational efforts. In practice, a trade-off is sought between accuracy and computational efficiency, which often renders the results of calculations qualitative, rather than quantitative. Furthermore, because the high dimensionality of the potential energy surfaces precludes a quantum mechanical description of the nuclear degrees of freedom, classical mechanics is used instead. Therefore, nuclear vibrations are not quantized and population transfers between the adiabatic states have to be modelled in an ad hoc manner with surface hopping, 7 or Ehrenfest dynamics. 8 In addition, the sampling of all relevant configurations may require more computational resources (and time) than is available. Nevertheless, the main advantage of atomistic models is that despite their limited accuracy, large and complex systems can be modeled directly, and provide qualitative insights into the effect of the chemical structure on the dynamics and energetics of a process.
The third term in Equation 1 describes the light-matter interaction within the dipolar approximation through g j (k z ): where u cav is the unit vector indicating the direction of the electric component of cavity vacuum field, here along the y-direction ( Figure S1); 0 the vacuum permittivity; and V cav the volume into which the mode with in-plane moment k z is confined. Figure S1: One-dimensional (1D) Fabry-Pérot cavity model. Two reflecting mirrors are located at − 1 2 x and 1 2 x, confining light modes along this direction, while free propagation along the z direction is possible for plane waves with in-plane momentum k z and energy ω(k z ). The vacuum field vector (red) points along the y-axis, reaching a maximum amplitude at x = 0 where the N molecules (magenta ellipses) are placed, distributed along the z-axis at positions z j with 1 ≤ j ≤ N .

One-dimensional periodic cavity
As in Michetti and La Rocca, 6 we impose periodic boundary conditions in the z-direction, and thus restrict the wave vectors, k z , to discrete values: k z = 2πn/L z with n ∈ Z and L z the length of the 1D cavity. With these approximations the molecular Tavis-Cummings Hamiltonian in Equation 1 can be represented as a (N + n mode ) by (N + n mode ) matrix with four blocks: The upper left block, H mol , is a N × N matrix containing the single-photon excitations of the molecules. Because we neglect direct excitonic interactions between the molecules this block is diagonal, with elements labeled by the molecule indices j: for 1 ≤ j ≤ N . Each matrix element of H mol thus represents the potential energy of a molecule (j) in the electronic excited state |S j 1 (R j ) while all other molecules, i = j, are in the electronic The |0 in Equation 5 indicates that the single-photon Fock states of all cavity modes are empty.
The lower right block, H cav , is a n mode × n mode matrix (with n mode = n max − n min + 1) containing the single-photon excitations of the cavity modes, and is also diagonal: for n min ≤ b ≤ n max . Here |1 b is the single-photon Fock state of cavity mode b with wave-vector k z = 2πb/L z . In these matrix elements, all molecules are in the electronic ground state (S 0 ), while cavity mode b is excited. The energy is therefore the sum of the cavity energy at k z , and the molecular ground state energies: Here ω cav (k z ) is the cavity dispersion (dashed curve in Figure 1b, main text): where ω 0 is energy at k z = 0, and c the speed of light.
The two N × n mode off-diagonal blocks H int and H int † in the multi-mode Tavis-Cummings Hamiltonian (Equation 4) model the interactions between the molecules and the cavity modes.
These matrix elements are approximated as the overlap between the transition dipole moment of molecule j and the electric field of the cavity mode b at the center z j of that molecule: for 1 ≤ j ≤ N and n min ≤ b ≤ n max .

Semi-classical molecular dynamics
The multiple cavity-mode Tavis Because the energies of the electronic ground and excited state, as well as the transition dipole where c m (t) are the time-dependent expansion coefficients of the time-independent polaritonic basis functions ψ m defined in Equation 11. A unitary propagator in the local diabatic basis is used to integrate these coefficients, 2,9 while the nuclear degrees of freedom of the N molecules evolve on the mean-field potential energy surface: Because the simulations in this work were run for only a few femtoseconds, cavity decay was neglected. Including decay into this model is, however, possible, as we have shown in previous works. 2,3 2 Simulation Details

Tetracene model
Prior to the MD simulations, the ground-state (S 0 ) geometry of the tetracene molecule was optimized at the CAM-B3LYP/6-31G(d) level of density functional theory (DFT), 10-13 while the excited state (S 1 ) geometry of tetracene was optimized at the TDA-CAM-B3LYP/6-31G(d) level. 14 A 10 ps MD simulation was performed for an isolated tetracene in the ground state at 300 K, using a stochastic thermostat with a relaxation time of 0.1 ps. 15 The equilibration simulation was used to compute the absorption spectrum shown in Figure 1b of the main text (see below).

Molecular dynamics of Tetracene-cavity systems
For the simulations in the single-mode cavities, the cavity mode was resonant with the tetracene excitation energy, which is 3.22 eV at the TDA-CAM-B3LYP/6-31G(d) level of theory. For simulations in the multi-mode cavities, the cavity was red-detuned, with a cavity energy of ω 0 = 2.69 eV at k z = 0 and cavity length of L z = 15 µm. The dispersion of this cavity was modelled with 60 modes (n max = 59), corresponding to an energy cut-off at 5.84 eV, which is sufficiently high above the S 1 energy of tetracene. Unless stated otherwise, the molecules were oriented to maximize the coupling strength by aligning their transition dipole moments to the polarization of the vacuum field inside the cavity at the start of the simulations (i.e., u y in Equation 10, Figure S1). In the multi-mode simulations, the molecules were distributed evenly along the z-axis of the cavity.

Rhodamine model
Prior to the MD simulations, we used density functional theory (DFT), [10][11][12][13] to optimise the groundstate (S 0 ) geometry of Rhodamine (Rh), shown in FigureS2, with the CAM-B3LYP functional 11,12 and the 6-31G(d) basis set. 13 The excited state (S 1 ) was optimized using time-dependent density functional theory (TDDFT), 20 within the Tamm-Dancoff approximation (TDA). 14 At this level of theory, the vertical excitation energy of Rh is 3.49 eV in the S 0 minimum geometry, while the energy gap to the ground state is 3.40 eV in the S 1 minimum.  Figure S1).

Molecular dynamics of Rhodamine-cavity systems
In this set-up, a Rabi splitting energy of Ω Rabi = 446 meV was achieved. The integration time step was 0.1 fs for all cavity simulations, which were performed with Gromacs 4.5.3, 16 in which the multi-mode Tavis-Cummings model presented above, was implemented, using the QM/MM interface to TeraChem. 17,18

Simulation Analysis
The absorption spectrum of the molecules was computed from trajectories in the ground (S 0 ) state.
These spectra were composed of the S 1 -S 0 energy gaps (1,000 snapshots) by super-position of Gaussian functions: 21 where I abs mol (E) is the intensity of the molecular absorption as a function of excitation energy (E), s the number of snapshots included in the analysis, ∆E i the excitation energy in snapshot i and µ TDM i the transition dipole moment. A width of σ = 0.05 eV was chosen for the convolution.
Following Lidzey and coworkers, 22 we define the "visibility", I m , of polaritonic state ψ m as the total incoherent photonic contribution to that state (i.e., I m ∝ nmax b |α m b | 2 ). Thus, the k z -dependent, one-photon absorption spectra of the strongly coupled molecule-cavity systems were computed as follows: For each frame of the trajectory, the polaritonic states were computed and the energy gaps of these states with respect to the overall ground state (i.e., E 0 , with all molecules in S 0 , no photon in the cavity: summed up into a superposition of Gaussian functions: Here, I abs (E, b) is the absorption intensity as a function of excitation energy E and in-plane momen- The photonic contribution of the total polaritonic wave function was computed as the sum of the projections of the cavity modes basis states (i.e., where N + n mode is the number of polaritonic states.  (Table S1). The evolution of the potential energy of the Tc molecules is plotted in Figure S3 for a bare molecule (i.e., outside the cavity) as well as for 1, 2, 4, 16, 32, 64 and 128 molecules inside the here to avoid obscuring the interpretation, the molecular dynamics are under damped. Therefore, the Tc molecule does not get trapped, but rather oscillates around the minimum on the S 1 potential energy surface instead. Nevertheless, the results of the simulations suggest that inside the cavity, the local S 1 minimum is accessible to molecules that are non-resonantly excited, irrespective of the number of molecules, or cavity field strength, in line with results from previous simulations. 27,28 Assuming, therefore, that dissipation will eventually dampen the dynamics and trap the excited Tc molecule in the S 1 minimum geometry, while the other molecules remain in their S 0 minimum geometries, we inspected the eigenstates of the Tavis-Cummings Hamiltonian in this combination of molecular configurations. The contribution of the molecule in the S 1 minimum energy geometry (molecule 1) to the lowest energy dark state, |β DS 0 1 | 2 (Equation 11), listed in Table 1, rapidly converges to unity upon increasing the number of molecules.

Analysis of non-adiabatic coupling vector
The non-adiabatic coupling vector connecting displacements in the vibrational modes of molecule i to transitions between polaritonic states ψ l and ψ m can be expressed as a Hellmann-Feynmann force: with ∇ a∈i the gradient with respect to the displacement of an atom a in molecule i and E l the adiabatic energy of polaritonic state l. After substitution of the expression for the polaritonic eigenstates ψ m (Equation 11), the Hellman-Feynman term becomes: 3 In what follows, we assume a single-mode cavity, i.e. n mode = 1. We also distinguish between dark states, which lack a photonic component (i.e., α m b = 0) and bright states, which contain a  Figure S4. With these assumptions, the Hellmann-Feynman term between the relaxed dark state (DS 0 ) and the LP, associated with a displacement of an atom in any of the molecules that are in the S 0 minimum, reduces to: In contrast, for a displacement of an atom in the molecule that is in the S 1 minimum (i.e., molecule 1) the Hellmann-Feynman term becomes: where the factor 0.9 was obtained from the fit to β LP 1 in Figure S4. This analysis suggests that the non-adiabatic coupling vector for transitions between the DS 0 and LP states is dominated by atomic displacements of the molecule that is in the S 1 geometry.
To verify this result, we performed MD simulations of 16 tetracene molecules, strongly coupled to a single-mode cavity with a cavity resonance ω cav = 3.22 eV and a cavity volume V cav = 15.6 nm 3 . The simulations were started in the DS 0 state with one of the tetracene molecules in the S 1 geometry, and all others in the S 0 geometry. Initial velocities were randomly selected from a Maxwell-Boltzmann distribution at 300 K. Two simulations were performed: in first simulation, the atoms of the molecule that is in the S 1 minimum (molecule 1), were kept fixed, whereas in the second simulation the atoms of all molecules that are in the S 0 geometry, were frozen. In Figure S5 we plot the population of the LP during the simulations. As shown in Figure S5a, there is no population transfer if the molecule that is in the S 1 geometry, is frozen. In contrast, when the other molecules are frozen, there is population transfer from the DS 0 state into the LP ( Figure S5b).
These results thus confirm that only displacements of atoms in the molecule that is in the S 1 state, can drive non-adiabatic transitions between the DS 0 and LP states. Figure S5: Population of the lower polariton (LP) in simulations of 16 tetracene molecules strongly coupled to a single-mode cavity, starting in the lowest energy dark state (DS 0 ) with one molecule in the S 1 minimum energy geometry, and all other molecules in the S 0 minimum energy geometry. In panel a, the molecule in the S 1 minimum energy geometry is kept frozen, whereas in panel b, all molecules that are in the S 0 minimum energy geometry, are frozen.
Because molecule 1 is in its S 1 minimum energy geometry, while all other molecules are in their S 0 minimum energy geometries (i.e., ∇ a∈1 V mol S 1 (R 1 ) = 0 and ∇ a∈j V mol S 0 (R j )=0 for j = 1), the

Hellmann-Feynman term in Equation 20
can be simplified further: Substituting this expression into equation 17 yields the non-adiabatic coupling vector connecting the lowest energy dark state to the lower polariton: The first term on the right hand side suggest that in order to mediate non-adiabatic population transfer between the DS 0 and LP, the molecular vibrations must have a geometric displacement, or Huang-Rhys factor, between the electronic ground (S 0 ) and excited (S 1 ) states. Because such modes also are Franck-Condon active, they show up in vibronically resolved absorption and emission spectra. The first term scales inversely with the square root of N . In contrast, the second term involving the cavity vacuum field and transition dipole moment of the molecules, is independent of N . Even if the number of molecules that fit inside the mode volume of a real Fabry-Pérot is finite, N can be very large and it is therefore unclear which of these two terms would dominate the non-adiabatic coupling vector in experiment. Nevertheless, as shown in Figure S6, both terms overlap with the same vibrational modes of tetracene. Therefore, to reach the main conclusion of this work that specific vibrations control the non-adiabatic relaxation from the dark states into the lower polariton, it is irrelevant which term dominates the non-adiabatic coupling. We note that because the energy gap between DS 0 and LP is proportional to the Rabi splitting, the prefactor also scales as 1/ √ N . However, for comparing different ensemble sizes in this work, we kept the Rabi splitting the same by scaling the mode volume of the cavity V cav with N : cav the mode volume used for simulations with 16 molecules in the cavity. Therefore, By keeping the Rabi splitting the same, also the energy gap remains the same. Hence, at constant

Relaxation in a strongly coupled Rhodamine-cavity system
In Figure Figure S8b shows this photonic weight (i.e., | m j c m α m j | 2 ) at 1 fs as a function of the vibrational energy of the mode along which the initial velocities were directed. The observation that population transfer predominantly occurs if we activate vibrational modes that overlap with the non-adiabatic coupling vector, suggests that also for Rhodamine, which lacks internal symmetry, the relaxation from the dark state manifold into the LP is selectively mediated by these vibrations.
As for Tetracene, we also computed the IR, Raman and vibronic spectra ( Figure S8). In contrast to Tetracene, not all modes that can drive population transfer, are Raman-active. Therefore, for Rhodamine, which lacks internal symmetry, the Raman spectrum is less useful for predicting such modes. Instead, these modes can only be identified from their overlap with the non-adiabatic coupling vector.

Model system for vibronic transitions between polaritons
To understand if in addition to overlap with the non-adiabatic coupling vector, the vibrational mode energy also plays a role in the relaxation from the DS 0 state into the LP, we constructed a simple model system. The model, illustrated in Figure S9, contains two polaritonic states, described by harmonic potentials. The minima of these potentials are separated in energy by ∆E LP,DS 0 , the gap between the LP and DS 0 states. Within the Born representation, the wave function of this model system is: where ψ m are the eigenstates of the Tavis-Cummings Hamiltonian (Ĥ TC ), here restricted to ψ LP and ψ DS 0 . The associated potential energy surfaces (V m ) are modeled as harmonic potentials: with where we omitted the ψ n |∇ 2 N |ψ m term, because it is much smaller than ψ n |∇ N |ψ m . The quantum number n l indicates the eigenstates of the harmonic oscillator in polaritonic state ψ l . The non-adiabatic coupling between the polaritonic states is conveniently evaluated using the stepping operators: whereâ|n = √ n|n − 1 andâ † |n = √ n + 1|n + 1 and we introduced the short-hand notation |n for |χ n . while in the second simulation, the energy gaps differed by 1 meV ( ω LP = E DS 0 −E LP ). Vibrational frequency and mass were obtained from the normal mode analysis, while the magnitude of the non-adiabatic coupling vector was extracted from the simulation with 512 tetracene molecules strongly coupled to the single-mode cavity. Figure S9b shows how the population of the vibrational states |0 DS 0 and |1 LP evolve during the simulations. Efficient population transfer from the ground vibrational state in the lowest energy dark state (DS 0 ) into the first vibrational excited state in the LP occurs efficiently only if the energy gaps match (top panel in Figure S9b). In contrast, when the gaps are different, transfer is inefficient (bottom panel in Figure S9b)